#include <stdio.h>
#include <sys/time.h>
#include <pthread.h>
#include <stdlib.h>
#include <math.h>
#include <errno.h>
#include <time.h>
#include <getopt.h>

static int SIZE;
static int NUM_THREAD;

struct option input_option[] = {
	{"size", 1, 0, 's'},
	{"thread_count", 1, 0, 't'},
	{0, 0, 0, 0}
};

struct init_param {
	int nIsTimeChecker;
	int nStartRow;
	int nLen;
	struct drand48_data gBuffer;
};
struct init_param* initParamList;

double **gMatrixA;
double **gMatrixB;

struct timeval startTime, endTime;

pthread_t* gThreadIDList;

void showMatrix();

pthread_barrier_t bar;

void* solveLinear(void* pData) {
	int i, j, k;


	struct init_param* pParam = (struct init_param*)pData;
	int start = pParam->nStartRow;
	int end = pParam->nStartRow + pParam->nLen;
	double dCoeffi;
	pthread_barrier_wait(&bar);
	if (pParam->nIsTimeChecker) {
	    printf("Start problem solving!!\n");
		gettimeofday(&startTime, NULL);
	}
	for(i=0; i<SIZE; i++) {
		for(j=start; j<end; ++j) {
			if (i == j)
				continue;
			dCoeffi = gMatrixA[j][i]/gMatrixA[i][i];
			for(k=0; k<SIZE+1; ++k) {
				gMatrixA[j][k] = gMatrixA[j][k] - dCoeffi*gMatrixA[i][k];
			}
		}
		pthread_barrier_wait(&bar);
	}
	for(j=start; j<end; ++j) {
		gMatrixA[j][SIZE] /= gMatrixA[j][j];
		gMatrixA[j][j] = 1.0;
	}
	pthread_exit(NULL); 
}

/**
 * Initialize program 
 * - Generate matrix A & array B
 */
int init()
{
	int i, j;

	/* Memory allocation */
	gMatrixA = (double**)malloc(sizeof(double *)*SIZE);
	for(i=0; i<SIZE; ++i)
		gMatrixA[i] = (double*)malloc(sizeof(double)*SIZE+1);

	gMatrixB = (double**)malloc(sizeof(double *)*SIZE);
	for(i=0; i<SIZE; ++i)
		gMatrixB[i] = (double*)malloc(sizeof(double)*SIZE+1);

	srand48(time(0));
	for(i=0; i<SIZE; ++i) {
		for(j=0;j<SIZE+1; ++j)
			gMatrixA[i][j] = gMatrixB[i][j] = drand48();
	}

	initParamList = (struct init_param*)malloc(sizeof(struct init_param)*NUM_THREAD);
	gThreadIDList = (pthread_t*)malloc(sizeof(pthread_t)*NUM_THREAD);

	return 0;
}

/** 
 * Finilize this program
 */
int fini()
{
	int i;
	for(i=0; i<SIZE; ++i) {
		free(gMatrixA[i]);
		free(gMatrixB[i]);
	}
	free(gMatrixA);
	free(gMatrixB);

	free(initParamList);
	free(gThreadIDList);
	return 0;
}

int gaussian_elim()
{
	int i, nQuota, nRemainder;
	int nNextRow = 0;

	nQuota = SIZE/NUM_THREAD;
	nRemainder = SIZE%NUM_THREAD;

	pthread_barrier_init(&bar, NULL, NUM_THREAD);
	for(i=0; i<NUM_THREAD; ++i) {
		if(!i)
			initParamList[i].nIsTimeChecker = 1;
		else
			initParamList[i].nIsTimeChecker = 0;
		initParamList[i].nStartRow = nNextRow;                                                                                                                                         
		initParamList[i].nLen = nQuota+(nRemainder-->0?1:0);                                                                                                                           
		nNextRow = initParamList[i].nStartRow + initParamList[i].nLen;
		if(pthread_create(&gThreadIDList[i], NULL, solveLinear, &initParamList[i])){
			printf("Error!!\t Fail to pthread_create()\n");
		}
	}
	for(i=0; i<NUM_THREAD; ++i){
		if(pthread_join(gThreadIDList[i], NULL)){
			printf("Error!!\t Fail to pthread_join()\n");
		}
	}
	printf("Finish problem solving!\n");
	gettimeofday(&endTime, NULL);
	return 0;
}

int getL2Norm(double* pL2Norm)
{
	double dL2Norm=0;
	double dTmp;
	int i, j;
	/* TODO Need to parallelize */
	for(i=0; i<SIZE; ++i){
		dTmp = 0;
		for(j=0; j<SIZE; ++j){
			dTmp += gMatrixB[i][j]*gMatrixA[j][SIZE];
		}
		dTmp -= gMatrixB[i][SIZE];
		dL2Norm += pow(dTmp, 2.0);
	}	
	*pL2Norm = sqrt(dL2Norm);

	return 0;

}

void set_option(int argc, char** argv)
{
	int opt;
	if(argc!=3 && argc!=5 ){
		printf("Usage: %s [options]\n", argv[0]);
		printf("\t-s\t--size\t\t Size of Matrix\n");
		printf("\t-t\t--thread_count\t Number of Thread\n");
		exit(0);
	}

	do{
		opt = getopt_long(argc, argv,"s:t:", input_option, NULL);

		switch(opt)
		{
			case 's':
				SIZE = atoi(optarg);
				break;
			case 't':
				NUM_THREAD = atoi(optarg);
				break;
			case -1:
				break;
			default:
				abort();
		}
	}while(opt!=-1);
}

int main(int argc, char** argv)
{
	set_option(argc, argv);

	double dL2Norm=0.0;
	init();
	gaussian_elim();
	getL2Norm(&dL2Norm);

	printf("\nL2Norm: %f\n", dL2Norm);
	printf("Time Ellapsed: %f msec\n", (endTime.tv_sec - startTime.tv_sec)*1000 + (endTime.tv_usec - startTime.tv_usec)/1000.0);
	fini();
	return 0;
}
